Dynamical Mean Field Theory of the Bose-Hubbard Model 
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Quantum dynamics of the Bose-Hubbard Model is in- 
vestigated through a semiclassical hamiltonian picture pro- 
vided by the Time-Dependent Variational Principle method. 
The system is studied within a factorized slow/fast dynam- 
ics. The semiclassical requantization procedure allows one to 
account for the strong quantum nature of the system when 
t/U <C 1. The phase diagram is in good agreement with 
Quantum Monte Carlo results and third order strong cou- 
pling perturbative expansion. 

PACS N. 74.20.-z, 05.30.Fk 



In recent years outgrowing interest has been devoted 
to many-body systems which can be modelled in terms 
of bosonic degrees of freedom. Examples are granu- 
lar superconductors, short-length superconductors and 
Josephson junction arrays [0. The relevant physics of 
these systems is captured by the Bose-Hubbard Model 
(BHM) which represents a boson gas of identical charges 
hopping through a D dimensional lattice. The boson dy- 
namics is described by the second quantized Hamiltonian 



■a-'jUi ) , (!) 
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where the operators = a]ai count the number of 
bosons at site i, while the annihilation and creation op- 
erators (Xj, a] obey the canonical commutation relations 
[dj, Oj] — Sij. The parameters of Hamiltonian (|l|) U > 0, 
t, and /i correspond to the strength of the Coulomb on- 
site repulsion, the hopping amplitude, and the chemical 
potential, respectively. The rich structure of the BHM 
phase diagram has been investigated by a number of the- 
oretical methods, from mean field variational H and 
perturbative Q approaches, to Quantum Monte Carlo || 
technique. 

At T=0 and integer fillings, the system undergoes to 
a quantum phase transition, between a Mott Insulator 
(MI) and a Superfluid (SF) phase. At t = the filling is 
fixed to the integer n that minimizes the onsite contribu- 
tion of the Hamiltonian (|l|). In this limit, the existence of 
a finite energy fi = 2U , required to add one boson to the 
system, reflects the MI nature of this phase, character- 
ized by a vanishing compressibility (gapped particle-hole 
excitations). The MI survives (except for the degener- 
ation points with /i/U = 2n ) when t/U > 0, inside 
extended lobes in the t/U, fi/U plane attached to the 



intervals I(n) = (2(n- 1)U , 2nU ") of the fi/U axis [§. 
Elsewhere, in the phase plane, the system is superfluid 
and appears to be compressible since the addition of a 
particle free to hop through the lattice brings up a re- 
duction of the Coulomb repulsion. A crucial property 
has to be mentioned, at this point, concerning the MI- 
SF phase lobe boundary. At the transition points the 
appearence of superfluidity is announced by the vanish- 
ing of the energy gap between the states corresponding 
to n (n — 1) and n+1 (n) particles (holes). Such a feature 
characterizes the whole frontier between the MI and the 
SF phase and will play a central role in the sequel. 

In spite of the great amount of work devoted to study 
the several aspects of the BHM, no investigation, to the 
best of our knowledge, has been made on the dynamics of 
boson degrees of freedoms. The attempt of relating the 
onset of the macroscopic order leading to the MI and SF 
phases with the microscopic behavior of the system mo- 
tivates the dynamical approach to the BHM. This should 
be useful as well to study both the emergence of vortex 
dynamics and the transport properties. In particular, in 
the present paper we shall show that the mean filed phase 
diagram of the BHM results to be quite improved if the 
dynamics of the SF order parameter is accounted for. 

Indeed it is possible to develop such a program, al- 
though several features of Hamiltonian (Q) (nonlinearity, 
the many-body character and its quantum nature) make 
the investigation of its dynamical behavior a hard prob- 
lem. 

The Time-Dependent Variational Principle (TDVP) 
method || offers a quite general procedure for construct- 
ing an approximate macroscopic wave function for many- 
body systems. Such a method, recently employed for 
studying the fermionic Hubbard Model dynamics @, is 
based on the idea of constraining the time evolution 
of the system's state |$) via the weaker form of the 
Schrodinger equation (Q\(ihd T — H)\&) = 0. Upon set- 
ting |$) = exp(iS /h)\Z) one obtains 



S = ih(Z\d T \Z) - (Z\H\Z) 



(2) 



where the trial macroscopic state \Z), the basic ingredi- 
ent of the method, must be structured so as to contain as 
much information as possible on the microscopic dynam- 
ics. The label Z is thus identified with a vector of micro- 
scopic parameters that are required to account for the 
dominating physical processes at microscopic level and 
represent the effective dynamical variables of the system. 
For this reason S and (Z\H\Z) can be properly seen as 
functions of Z and we denote them by C[Z] and H[Z], 
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respectively. The TDVP theory identifies the effective 
Lagrangian of the system with C[Z] and its semiclassical 
model Hamiltonian with TC[Z]. The quantum dynamics 
is thus described through a set of semiclassical hamilto- 
nian equations obtained by implementing the stationarity 
condition 5S = on the action S of C\Z\. 

A thorough analysis of the semiclassical equations re- 
lated to H is interesting in itself due to their complexity. 
Here, we intend instead employ the TDVP method for 
finding the ground-state configurations of the model and 
examining their dependence from the microscopic param- 
eters t/U and n/U. In other words, we aim to work out 
the phase diagram of the BHM starting from the study 
of the semiclassical equations of motion. Let us apply 
the TDVP method to the BHM model (§J). We write the 
state |$) as |$) = exp (iS/h)®i \zi) once the trial macro- 
scopic state \Z) is assumed to have the form \Z) = 
Here the states \zi) are the Glauber coherent states || as- 
sociated with the boson lowering operators aj, solutions 
of the equation a,i\zi) — z.i\Zi) for each i. The choice of 
\Z), of course, is suggested by the fact that H belongs to 
the enveloping algebra of the A s -boson Weyl-Heisenberg 
algebra {I, aj,a,j,rij : j G A}, N s being the number of 
sites of the lattice A. In this case Eq. (||) becomes 

C[Z] = -(z~iZi — ZiZi) - H(Z) , (3) 

i 

where the semiclassical model Hamiltonian H(Z) 
XII /. is easily shown to have the form 

n = ]T(t/N 2 - ^)M 2 - \ Yu + • ( 4 ) 

• (id) 

Lagrangian (|j|) yields the equations of motion 

ih£i = -/j,Zi + 2Uzi\zi\ 2 — - ^ Zj , (5) 

ieCO 

where (i) indicates the set of the nearest neighbour sites 
around i and we have omitted the equation for z% di- 
rectly ensuing from Eqs. (^) via complex conjugation. 
Notice that Eqs. (||) can be obtained as well through the 
standard formulas ifiZj = {zj,7i}, based on the canoni- 
cal Poisson Brackets {zk, Zj} = Skj/ih replacing commu- 
tators [cij,<xt] = Sij in the TDVP semiclassical scenery. 
This checks their hamiltonian character. As expected, 
Eqs. (j^) are not integrablc (integrability occurs only 
when either U — or when t = 0) since the only known 
constant of motion, a part from TL, is the semiclassical 
version M = J^- \zi\ 2 of number operator N = J^. rij. 

We simplify the structure of Eqs. (||) by separation of 
slow and fast dynamics, a procedure which is in some 
way the analog, in a dynamical contest, of the Mean 
Field Approximation (MFA) usually employed in Sta- 
tistical Mechanics. To this end we set, at each site, 
Zi = tpi + rji and assume that ipi is a slow variable 



whereas r\i is a fast oscillating term describing the com- 
plex, high-frequency part of the dynamics taking place 
on the microscopic interactions time scale (the hopping 
interaction, in this case). Then ipj — (zj) T ((•) r de- 
notes time average) when the time scale r is larger than 
that of the rjj's. The onset of order in the system at 
the macroscopic scale should reflect the dominating role 
of the 'ipj's in the lattice dynamics. Imposing the stan- 
dard condition (zi — ipi){zj — 'ipj) = rjifjj w of the MFA 
procedure involves the dynamical mean field decoupling 
ZiZj ps ipiZj + tpjZi — ipiipj which implies, in turn, the 
random phase approximation (zjZi) T w (zj) T (zj) r . The 
dynamical scenery just depicted together with an ergodic 
assumption leads thus naturally to defining 

s j 

as the macroscopic order paramenter revealing when or- 
der issues from the lattice dynamics. 

When our dynamical MFA is applied to Hamiltonian 
(^), and the further assumption ipj = ipi for j G (i) is 
made (smoothing condition), the kinetic term modifies 
as follows 

(id) i 

where q denotes the number of nearest neighbours per 
site. The resulting Hamiltonian reduces thus to a sum of 
on-site terms: H m f = J2j where 

Hj = U\z 3 \ 4 - n\zj\ 2 - | (zjiPj +iijZj - |^f) • (8) 
The hamiltonian equations ensuing from Eq. (Q) 

ihii = -[iz t + 2Uz i \z l \ 2 - y ipi , (9) 

bear memory of the off-site dynamics only through the 
on-site term ipi . When compared with the exact ones, (|^), 
they imply the relation qipi w Yljett) z j consistently lead- 
ing to an identity once time average is carried on and the 
smoothing conditions are used. It is important to notice 
that Af — J2i | | 2 is no longer a constant of motion while 
the mean field Hamiltonian 7i. m f = Ylj > that should 
approximate TL, might represents a time-dependent total 
energy (which Ti. is not) since it depends on i/^ 's. In order 
to recover such features, it seems reasonable to introduce 
some restriction on the form of ipj(r), whose time be- 
havior, so far, has not been specified at all. We do this 
looking for solutions of Eqs. (g) where 6j, \j, the phases 
of Zj = \zj\e l9: > , ipj = \ipj\e lXj respectively, are locked one 
to the other in such a way that A(p == Oj — Xj = const, the 
constant being zero or ir. In view of Eqs. (0), this entails 
that any \zj | 2 has a zero time derivative thereby restoring 
the particle total number M to its proper role of time- 
independent quantity. The same effect is obtained for any 
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Tij , and thus for Tt m f , only adding the further condition 
\tpj\ = const. Solutions of the form Zj(r) = \zj\e l6j ^ 
whose phase, due to Eq. (0), obeys the equation 



h\zj 



(2U\z 3 \ 2 



u)\ Zj 



qt 



(10) 



successfully fulfill the conditions just stated. Despite the 
elimination of any residual dynamical complexity we are 
able to characterize the MI and some features of the SF 
phase via the phase dynamics of Eq. (|lC|) . 

We examine first the dynamics related to the MI. In 
this case, ipi must have a zero time average along macro- 
scopic time scales. Such a behavior occurs when the 
uniform filling conditions rii = n, for all i (we identify 
here number operators n^s with their integer spectral 
values) is inserted in ( |To| ) by setting \zi\ 2 = n. Such a 
substitution is the natural consequence of the requanti- 
zation process Q of the action-like variables \zi\ 2 (notice 
that {|zi| 2 , 9j} = 8ij/h) strongly re ques ted from the pure 
quantum character of the MI. Eq. (|10|) is thus solved by 
0j(r) — X±r/h + ctj, (ctj is the initial condition) where 



\/n \ 2 



(11) 



S = u/U — 2n , and A_ (A + ) is related to the choice 
A<p = 7r (A<fi = 0). Notice that the index j does not label 
A± since the request (zj) T = ipj(r) leads to = y/n at 
each site. In the present theory, the frequencies X± play 
the role of time correlation lenght governing the phase 
transition. Our theory gives A± = U^/n(p — u c ) for fixed 
t and A± = q\ipi\/2{t — t c ) for fixed \i (u c and t c are the 
critical values of u and t) . Defining the critical exponents 
z and v as in the ref. S, we argue that M 



1. 



(12) 



By replacing in the reduced Hamiltonian (||) the value of 
provided by Eq. (fll), the energy of the MI reads 



£„(",£; A±) = n 



U5 -Un+ — (X± + U5) 2 
it 



(13) 



where the subscript n reminds us that the filling n 
is accounted for. The oscillating behavior of ^ = 
(e lX±T /N s ) J2j ipje iaj , having a vanishing long time av- 
erage, identifies the MI. This, in fact, implies that the 
gauge symmetry breaking expected in the SF phase can- 
not take place. Notice that the ordinary (time indepen- 
dent) MFA cannot describe the MI for t > 0, since the 
hopping term of the reduced Hamiltonian is canceled by 
the vanishing of the order parameter, -0 = 0. Within our 
scheme, instead, the condition (^) T — can be realized 
also for ip ^= 0. 

All dynamics disappears at the degeneration points 
u/U = 2n, t/U = 0, that represent the limiting points 
of the SF domain separating I(n) from I{n + 1). Within 
our theory, the fixed points of the equations of motion 



(A± = 0) correspond to the SF system configurations. 
This is but the oversimplified version of the low frequency 
dynamics expected in the SF phase. Considering such 
configurations (the trivial CclSG Zj = due to Zj — ipj = 
is excluded) allows one to recast Eq. (||) in the form 



(14) 



making tpj a function of Zj . The energy associated to the 
hamiltonian Tij , in turn, reduces to 



e((i, t, Zj 



M 2 



qt 



(u-2U\z J \ 2 y + (u-3U\z J \ 2 ) 



once ( |i4| ) is inserted in (||). e(u,t,Zj) is the on-site en- 
ergy accounting for the absence of dynamics. The limit 
X± — > 0, in fact, shows that E n (u,t; X±) — > e(fx,t,Zj) 
provided \zj\ 2 — n. Its lowest value is found to be, via 
minimization, e(£) = — J7|C| 4 , where the value of \zj\ 2 cor- 
responding to the minimum is |£| 2 — (u + 2t) / (2U) . It is 
worth noticing that inserting |£| in ( |14| ) implies tpj = Zj 
so that the minimum energy configuration naturally sat- 
isfies the consistency condition on which our dynamical 
MFA is based. 

Let us consider, instead, how the MI characters reflect 
on the macroscopic phase - viable to experimental obser- 
vations - S. Notice, first, that inserting Eqs. (||) in the 
Lagrangian (||) implies S = U J2j \zj\ 4 and that the same 
result is found in the MFA scheme in that ipj « Zj . Then 
the exponential factor of |<I>) has a frequency S which 
reduces to 



S = UN s n 2 



(15) 



when the system is a MI with \zj\ 2 = n for each j. A 
transition leading from the n-lobe to the (n+l)-lobe thus 
involves a change of the phase frequency of UN s (2n+ 1), 
whereas in case of a transition between any two superfluid 
states, where S m f — UN S \^\ 4 , since Zj = ^ for each j, no 
quantization of the frequency variation appears, in that 
^ takes continuous values. 

Now, we employ the expression ([l3]) for the on-site en- 
ergy to reconstruct the lobe-like structure of the phase 
diagram. In the SF phase, the states with n and n + 1 
(adding a particle) , as well as the states with n — 1 and 
n (adding a hole) must be degenerate. The curves repre- 
senting the n-lobe boundary are identified by implement- 
ing both gauge symmetry breaking through the limits 
A± — > and vanishing of the energy gaps E n — E n ±\. In 
other words we require 



lim (E n 



E n+1 ) = (6 < 0) 



lim (£?„_!- E n ) = 



(6>0) 



(16) 



(17) 



For solving Eqs. ©-(0) 

we introduce the variables 
5± = u/U - 2n + (1 ± 1). By inserting 6+ > (<5_ < 0) 
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in eq. (|T^) ((0)), and defining r — qt/4U, one gets the 
quadratic equations 5± + 2r8± — 2r (2n=pl) = providing 
the curves 

^ = 2(n - (1 ± 1)) - r ± [r 2 + 2r (2n T 1)] 1/2 • (18) 

The lower branch and the upper /J,-(t) constitute 

the boundary encircling the n-th lobe. We conclude by 
retrieving from Eq. (jl^) the position of the farthest 
point on the n-lobe boundary from the /it-axis. By setting 
H-{t) = one finds the lobe tip coordinates 

t = U/qn (19) 

and n(t c )/U — 2n — 1 — (l/2n). In the captions of Figs. 
(1) and (2) the values of t c furnished by the present ap- 
proach, is compared for the D = 1 and D = 2 cases with 
QMC H and the Strong Coupling Perturbative Expan- 
sion (SCPE) @. 
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FIG. 1. The Phase Diagram of the BHM for D = 1. The 
error boxes are the QMC results of Batrouni et al. Rel- 
atively to the first lobe (n, = 1), (t c /U) = 0.5. QMC 
gives (t c /U) = 0.43 ± 0.002 and the SCPE (t c /U) = 0.43. 
For (n< = 3), QMC and SCPE give (t c /U) = 0.2 and 
(tc/U) = 0.18 respectively. Our theory gives (t c /U) = 0.16. 

The dynamical approach we have developed appears 
to succeed in describing the quantum MI-SF phase tran- 
sition of BHM. The resulting phase diagram indeed ex- 
hibits an excellent agreement with QMC simulations and 
SCPE results. This suggests that Eqs. (||), here faced 
solely within the dynamical MFA, deserve a systematic 
investigation by the methods of dynamical system theory. 
The dynamics they account for, in fact, should describe 
not only zero temperature configurations but also excited 
states involving density waves as well as vortices. More- 
over, a systematic analysis of Eqs. (0) should be interest- 
ing both in relation to the dynamical scaling theory || 
and for calculating the dynamical correlation functions 
in the MI. Work is in progress along these lines. 
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FIG. 2. The Phase Diagram of the BHM for D = 2. 
The error box indicates the QMC tricritical point obtained 
by Krauth and Trivedi. For m = 1, (t c /U) = 0.25 while 
QMC gives (tc/U) = 0.244 ± 0.002 and SCPE provides 
(tc/U) = 0.272. 
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